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Abstract — We present the sensor setup and the basic naviga- 
tion algorithm used for the flight control of the SkySaUs towing 
kite system. Starting with brief summaries on system setup and 
equations of motion of the tethered kite system, we subsequently 
give an overview of the sensor setup, present the navigation 
task and discuss challenges which have to be mastered. In 
the second part we introduce in detail the inertial navigation 
algorithm which has been used for operational flights for years. 
The functional capability of this algorithm is illustrated by 
experimental flight data. Finally we suggest a modification of 
the algorithms as further development step in order to overcome 
certain limitations. 

I. Introduction 

The SkySails system is a towing kite system which allows 
modern cargo ships to use the wind as source of power in 
order to save fuel and therefore to save costs and reduce 
emissions [1]. The SkySails company has been founded in 
2001 and as main business offers wind propulsion systems 
for ships. Starting the development with kites of 6-10 m^ 
size the latest generation of products, with a nominal size 
of 320 m^, can replace up to 2MW of the main engine's 
propulsion power. Besides marine applications there are 
strongly increasing activities in using tethered kites and rigid 
wings for generating power from high-altitude wind [2], [3]. 
Therefore the design of control systems for tethered kites 
has become a growing field of experimental [4], [5] and 
theoretical [6], [7], [8], [9] research efforts. These control 
designs demand stable and reliable state estimates under 
changing wind conditions in order to perform their complex 
task in a satisfying and robust way. 

In this paper we focus on the sensor setup and on the 
estimation algorithms used for our flight controller. The 
most important task is the estimation of the kite orientation. 
In contrast to many applications in the field of inertial 
navigation we do not have a sensor which provides absolute 
attitude information. The task of attitude estimation, based on 
acceleration and rotation rate sensors alone, poses a special 
challenge for our navigation algorithm. 

Many classical inertial navigation algorithms are based on 
Kalman filtering [10], [11], [12] and extensions to this theory 
[13]. Our approach is related to complimentary filtering [14], 
[15] which will be discussed in the last section. We would 
like to emphasize that the aim of this paper is to illustrate 
the principle ideas and the challenges to be tackled due 
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Fig. 1. Definition of coordinates for the considered system. The right- 
handed coordinate system is defined by the basis vectors e^ , e^ , e^ with 
e^ in wind direction and e^ pointing downwards with respect to gravity. The 
kite position is parameterized by introducing two different sets of spherical 
coordinates {(/>, 5} (left figure) and {ip, i9} (right figure). For a more precise 
definition see J2| and |3l respectively. The kite axes are labeled as roll 
Erolh pitch Gpitch and Y^w Byaw TMs Corresponds to the definition usually 
used in aerospace applications with roll axis parallel to forward and yaw 
axis parallel to down directions respectively. Note that for the usual flight 
situation the yaw vector eyaw is defined by the position of the kite assuming 
it is more or less constrained to the origin by a rigid rod. Thus orientation 
of the kite is represented by the single angle i/ig and i/i respectively. 



to system dynamics and environmental disturbances specific 
for our combination of sensors and towing kites. Thus the 
presentation also reflects the development history and does 
not intend to be a rigorous, ready-to-use construction manual 
for a general type of filter algorithm. 

The paper is organized as follows: We start with short 
summaries of the overall system setup, the selected coor- 
dinates and the system dynamics. We then describe our 
sensor setup and present the attitude estimation algorithm in 
detail. Subsequently we explain the wind referencing setup 
and justify the navigation function by experimental flight 
data. We complete the article with a discussion of further 
development steps extending the previous algorithms. 

II. System setup and Coordinates 

In the following we will make use of two different 
coordinate systems for the navigation and control algorithm, 
respectively. The control design is based on a coordinate 
system with the wind direction as symmetry axis, while 
the attitude estimation algorithm is referenced to the gravity 
(down) direction. The choice of these different coordinate 
systems reflect the difference in the two design tasks (see 
[16] for further discussion especially of the controller related 
aspects). 

The two coordinate systems are introduced in Fig. [T] For 
a line length L the state of the kite is defined by the set of 
the three angles {0, 6', ^/ig} and {(/3,i9,-0} respectively. With 
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Fig. 2. Principle of system dynamics: tlie kite position is given by i?, ip 
(left figure) wliile the yaw angle ip denotes the orientation with respect to 
the wind. The quantity ip is crucial for controlling the flight direction of 
the kite, shown for the figure-eight pattern (right figure). 



respect to the basis vectors e^, Cy, Cz the kite position x is 
determined by the set {(j), 9} and {(p, ^} as: 

cos cos \ / cost? 

X — L \ sin (^ cos 9 \ = L \ sin ip sin i? 

— cos ip sin d 



(1) 



— sm ( 



The kite axes are denoted as eroii (roll or longitudinal), cTpitch 
(pitch) and Cyaw (yaw). 

For a description based on rotation matrices one would 
start with a kite at position Le^, with roll-axis in negative z- 
direction eiou = —ez, and then apply the following rotations: 
first —ipg about x, then 9 about y and finally (p about z. This 
transformation reads: 



R^Rz{(j))Ry{9)R^{-^^ 



(2) 



The quantity %p^ represents the angle between the Cpitch-axis 
and the horizon or the angle between the eroii-axis and the 
'upward' (— e2)-direction. 

In the second coordinate system {ip, •&, ip} the rotations 
are first —ip about x, then t? about y and finally ip about x. 
This transformation reads: 



R^R,i^)Ry{^)R,i-ijj). 



(3) 



Here one could interpret the angle ip as orientation of the 
kite longitudinal axis with reference to the wind. For a given 
kite position x (parameterized by ip and i?) the reference 
orientation ^ = corresponds to the minimum of the scalar 
product (eroiijCa;) obtained when turning the kite fixated at 
this position x around its yaw axis eyaw A nonzero value of 
ijj represents a kite orientation obtained by a rotation of tp 
about the yaw axis Cyaw starting at this reference orientation. 
Next we provide the transformation between the two 
coordinate systems {(t),9,ipg} ^ {ip,d,ip}: 



iP 



arctan(sin (p cos 9, sin 9) 

arccos(cos(/)cos6') 

■^g + arctan(sin p, cos (p sin ( 



(4) 
(5) 
(6) 



III. Dynamics and Equations of Motion 

For the description of the dynamics we refer to Fig. l2] 
illustrating a dynamically flying kite in a figure-eight pattern. 
The position of the kite is parameterized by the previously 



TABLE I 

Overview of sensors aboard the control pod (upper rows) 

and the ship (lower rows). 



Sensor 


Quantity 


Description 


Inertial Measure- 


oJs 


Turn rates of the control pod below the 


ment Unit (IMU) 




kite 




Qs 


Accelerations of the control pod below 
the kite 


Impeller 


Va 


Air path speed in ejon -direction 


Anemometer 






Strain Gauge Pod 


F 


Force measurement towing line 


Barometer 


h 


Barometric height of the control pod 


Tow Point 


0s, 9b 


Angles of towing rope with respect to 
the ship (compare Fig. T left). 


Ship Anemome- 


fw. 't>vi 


Apparent wind speed and direction 


ter 




aboard the ship. 


Ship IMU 




Roll and pitch angles of the ship. 


Line length 


L 


Rotary encoder measurement of un- 
wound towing line length 



introduced coordinates ip, d while the flight direction is 
determined by -p. Before explaining the dynamics we give a 
summary of the equations of motion: 



iP 



gVg^S + ip cos 'd 

Ua f , tan?9 
-(^cos^-^ 



Lsin^ 



sin?/). 



(7) 
(8) 

(9) 



The basic kite response is given mainly by the first term of 
(JTJi while ip cos-d is a correction due to the motion on the 
sphere [16]. A steering deflection S results in a turn rate about 
yaw axis, scaled by g and air path speed Ua. The integrated 
quantity ip represents the kite orientation with respect to the 
wind and thus controls the flight direction (compare (J8]l and 
Q). E denotes the glide ratio of the kite. For an angle 
-0 = the system stays more or less stationary with the 
kite roll axis aligned against the wind direction. In order 
to obtain the depicted dynamical 'figure eight' pattern one 
has to command a certain -08 > 0, resulting in a kite motion 
into negative ej,-direction. After some time the commanded 
i/jg is changed to a negative value leading to a curve and 
subsequently to a motion into positive ej,-direction. 

The altitude -d of the flight trajectory can be controlled 
by the amplitude of \ips\. We would like to stop at this 
point and not to proceed into a detailed derivation. We just 
conclude that ip is the crucial quantity for the flight controller. 
A detailed justification of the equations of motion and the 
control system architecture can be found in [16]. 

IV. Sensor Setup 

An overview of the sensor setup is given in Table |l] 
The data acquisition is performed by distributed computers 
running and merging the data at a main sample rate of 10 Hz. 

The main sensor is the inertial measurement unit (IMU) 
in the control pod, which is situated under the kite including 
three perpendicular arranged gyroscopes and accelerometers 



based on the solid state microelectromechanical systems 
(MEMS) technology. For some years we used the Crista 
OEM sensor head from Cloudcap Technology Inc. [17] 
which combines sensors from Analog Devices Inc. with 
an A/D-converter and includes a factory calibrated lookup 
table for gains and offsets over the operational temperature 
range. Later we switched to the ADIS 16355 device from 
Analog Devices Inc. [18] which provides a complete IMU 
with temperature calibrated outputs on the digital interface. 

We would like to emphasize at this point that calibration 
and offset compensation over the operational temperature 
range are crucial for the performance of the navigation 
algorithm. In order to improve the acceleration offsets of 
the ADIS devices we determine the offset values on a water 
leveled reference plate before installation of the sensor units. 
Both IMUs have a measuring range of ±300deg/s (gyros) 
and ±100 m/s^ (accelerometers). The typical observed bias 
values (Icr) for both IMUs are in accordance with the data 
sheets and approximately 0.2 deg/s for the gyros and 0.2 m/s^ 
for the accelerometers. 

The IMU values are sampled at a rate of 200 Hz. 20 values 
are averaged for the turn rates cDg and for the accelerations 
as which are then provided to the main 10 Hz cycle. 

The air path speed of the kite in the range of 5-50 m/s is 
determined by a frequency measurement of pulses generated 
by an impeller anemometer 

For completeness we note that there are also a strain gauge 
and a barometer altimeter aboard the control pod. They are 
used for redundancy and guidance purposes. However, these 
topics will not be elaborated in this paper 

The towpoint is the deflection point of the towing rope 
aboard the ship. Two angular sensors determine the mechan- 
ically sensed direction of the towing line with respect to the 
ship as azimuth and elevation angles. Another IMU device 
(XSens MTi) is installed on the ship measuring the roll and 
pitch angles due to the wave induced motion of the ship. 

The azimuth and elevation angles are merged with the 
roll and pitch angles in order to get the azimuth (p^ and 
elevation 9^ values. This corresponds in good approximation 
to an inertial reference and provides reasonable suppression 
of sea state influences which could easily cause ship roll 
amplitudes of ±10deg in operating conditions. We would 
like to point out that above transformation is correct only if 
the ships' center of angular motion coincides with the tow 
point, which is usually not the case. Yet, it is an acceptable 
approximation as the main dynamics takes place with respect 
to the roll axis. 

A further sensor is the anemometer aboard the ship 
measuring wind speed and direction of the apparent wind. 
Although we will argue later that the wind conditions at flight 
altitude might be quite different, this sensor is nonetheless 
an important reference input for the estimation algorithms 
especially during the launch of the system. 

The dynamics is influenced by the line length which is 
measured by an multiturn rotary encoder attached to the drum 
of the towing rope winch. 

For our prototype systems we also used two GPS receivers 
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Fig. 3. The yaw estimator's state is represented by the two rotations Rs 
and i?r ■ Rotation Rs is computed by quaternion integration of the three turn 
rates and depicts the fast dynamics of the flying system. The rotation iJr 
references the system to the inertial frame thus capturing a potential slow 
drift due to sensor errors. 



in the control pod and aboard the ship for calculating the 
difference between positions. Although we obtained valuable 
data for research and development purposes we discarded the 
use of GPS for operational flights as the position data suffer 
from dropouts due to rapid changes of antenna orientation 
during dynamic flights. 

V. Yaw Angle Estimator 

In this section we discuss the yaw angle estimator (YAE) 
which determines ^/jg. This quantity can be considered as 
the attitude angle relative to the horizon. It is important to 
emphasize that we have no direct measurement of this angle, 
e.g. based on optical or microwave technology. Such a system 
would be expensive and difficult to operate at sea. Therefore 
we must estimate the orientation based on the rotation rates 
Wg and the acceleration values a^ alone. 

The orientation can be determined by integration of the 
turn rates but has to be referenced to the horizon. In order 
to obtain the direction of gravity we carry out the following 
steps: first the measured accelerations are transformed to the 
integration reference frame and then the obtained accelera- 
tion vector is averaged. This coordinate system represents, on 
average, an inertial frame: As we deal with a tethered system 
we expect that due to the averaging all dynamic accelerations 
cancel out over time, because otherwise the system 'would 
fly away' in contradiction to the constraint by the tether 
Thus the average value represents an estimate for the gravity 
vector However, extra precautions have to be taken in order 
to deal with sensor offsets, which may lead to a slow drift 
of the coordinate system. 

The main idea of the algorithm is illustrated in Fig. l3] 
the turn rates ajg are integrated and then used to build an 
orientation matrix i?s, which represents the orientation of the 
IMU sensor with respect to the integration reference. This 
reference already features the property of an inertial system 
but drifts due to the measurement errors and sensor offsets 
with respect to the inertial reference frame. By evaluation of 
accelerations and extraction of the apparent gravity direction 
a rotation matrix R^ is computed by averaging. Thus R^ 
relates the, possibly drifting, integration system to the inertial 
reference frame. The final navigation angles are calculated 
on the basis R — R^-R^ which represents the orientation of 
the IMU sensor with respect to the inertial reference system. 
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Fig. 5. Setup of the navigation algoritlim wliich references the 0g output 
of the YAE and computes the orientation tpm with respect to the wind for 
the flight controller. 



Fig. 4. Overview of the yaw angle estimator algorithm which computes 
the orientation estimate {<j>g, 9g, '/'g} using turn rates tUg and accelerations 



We now will go into the details of the estimator algorithm. 
The respective components and their interactions are shown 
in Fig. HI The first block to be mentioned is the integration 
of the turn rate vector uj = LUg — uJo to Rg. The vector cjq 
denotes the sensor bias estimates. Details of the quaternion 
propagation can be found e.g. [19]. The rotation matrix Rg 
describes the fast dynamics of the systems. 

This matrix Rg is now used to transform the measured 
accelerations ag into the integration reference by 

(10) 



*avg 



This vector is subsequently filtered by a lowpass: 
LP[dr], where LP[] represents three independent filter op- 
erations on the vector components using a second order 
Butterworth lowpass with a cutoff frequency of 0.01 Hz. 



We would like to emphasize that ( 10 1 and the subsequent 
filter realize the central idea introduced in the beginning of 
this section: The component wise averaged value Savg thus 
represents an estimate for the gravity vector 

The next task is to compute the matrix R^ such that 
the normalized acceleration estimate Sq = a. 



avg/ 



-*avg| 



transformed into the (—e^) -direction which is the measured 
acceleration signal of the gravitational force pointing into 
'down' -direction: 

-RrflO = -4- (11) 

The straight forward solution of the R^ propagation task 
could be to determine the drift of this vector between time 
steps n — )■ (n + 1) by Aa — ao{n + 1) — ao{n) and 
propagate rotation R^ accordingly. Since numerical errors 
would accumulate a better solution is the introduction of a 
feedback loop and basing the propagation upon the deviation 

Aa — do ~ Rj. 



(-e.). 



(12) 



We then have to find a rotation fl which compensates for 
this deviation by turning Sq accordingly: 



Aa = —Q X clq. 



(13) 



Inversion of this equation with respect to fl is not unique. If 
we require fi _L ap, which is equivalent to the minimal rota- 
tion fulfilling (fTSll, some elementary vector algebra yields: 



This rotation f2 is used for the propagation of R^ via 
quaternion algebra. The navigation angles {0g,6'g,V'g} ^"s 
calculated out of i? = Ri-Rs by using a similar sequence of 
rotations as given in (pi. 

For perfect sensors we would have finished the task at 
this point, but our gyroscopes bring along bias values on the 
turn rates which would accumulate to angle errors. In order 
to minimize these angle errors, it is crucial to compensate 
for these offsets before the quaternion integration (compare 
Fig. |4|. Explaining the estimation of the offset vector lDq we 
start with the error due to dragging which is defined as: 

Aar = a,. - aavg- (15) 

Now a rotation (l^ has to be found with the property 

ASr = r^i- X a,.. (16) 

It is worth mentioning that this is a rather intuitive approach 
as during curve flights |ai-| ^ I'J^avgl- Therefore the relation 
to a rotation seems to be weak. But as (l^. is strongly filtered 
the approach holds and worked in practice. An important 



issue of ( 16 1 is that Oi is the reference for Aor and not aa 



in order to get the smoother behavior With the same type 
of reasoning as for ([T4]i we get 

n, = -—2 (a, X Ad,) . (17) 

I a,. I 

For estimation of the offset rate ujq the transformed rotation 
R^Cli is integrated with a gain of 7 = 0.003 1/s as follows: 



cDo(n + l) =a;o(") + At7i?f^r• 



(18) 



n 



-dg X Ad. 



(14) 



We are aware that this offset estimation algorithm is more 
an intuitive than mathematically rigorous approach. In our 
description we further omitted an additional damping term, 
regarding the rotational degree of freedom with respect to 
6*2, which we have also implemented in order to further 
robustify the filtering scheme. Explaining these technical 
details cannot be be done within the limits of this paper. 

VI. Wind Referencing 

As akeady introduced, the controller input value ^/im rep- 
resents the orientation with respect to the wind. The general 
arrangement for its computation is shown in Fig. [5] The IMU 
turn rates Wg and accelerations dg are fed into the yaw angle 
estimator (YAE) which determines estimates for the angles 
9g and V'g with respect to gravity. As this gravity reference 
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Fig. 6. Two different dynamic regimes have to be covered by the navigation 
algorithm. In the left figure the regular constrained flight situation is shown. 
The towing line force orients Cyaw in line direction. For the free flight 
situation in the right figure the orientation of the kite has to be considered 
independently from the position. 



is a single vector only, the 0g output represents the dynamics 
with respect to the vertical Cz-axis but lacks an absolute 
reference and therefore shows an offset drift due to sensor 
noise. As a consequence the 'reference block' is introduced 
in order to reference this value to <j>i ^ (ps + {t^ ~ 4>w), 
which is the horizontal angle of the towing line to the 
wind direction, calculated from towpoint angle 4>s and wind 
direction 0„ aboard the ship. Without going into details this 
block performs a complementary filtering similar to [14] with 
highpass behavior with respect to ip^ combined with an offset 
generated by lowpass filtering of 0^. 

The computation of tpni is based on (f>gr, 9g and ipg using 
(|6]l. The specific choice of these input quantities is explained 
with Fig. |6] illustrating the different flight situations. In the 
left figure the regular mode of operation is shown which 
is determined by a constrained dynamics. The eyaw-vector 
points into direction of the towing line. Both the IMU and the 
towpoint measurement values coincide 0gr~0s + (tt — 0w), 
6g sa 6s and thus can be used likewise in order to compute 
ipni for the control task. 

A more challenging situation is the free flight situation 
shown in the right part of Fig. |6] Sudden wind fall offs, 
downwashes or wave induced motion of the ship may lead 
to flight situations where the distance between tow point and 
kite is shorter than the towing line length. These situations 
lead to a completely different dynamics of the flying system 
which now more or less behaves like a paraglider. Although 
in most cases these situations fortunately last only a few 
seconds, arrangements must be made that the input values 
for the flight controller stay within appropriate range. 

The computation of tp^n from the data set {(pgi-,6g,ipg} 
is based on control pod sensors apart from the referencing 
of (pg,- to ship sensors. As this referencing takes place on a 
long-term timescale the short-term timescale navigation can 
be considered as a local navigation aboard the control pod 
which is presumably the best orientation estimate of the kite 
with respect to the wind even in extreme flight situations. 

We conclude this section by noting that in order to get 
high-performance navigation information the determination 
of the wind direction is a crucial challenge as wind direction 
at flight altitude might differ from the wind direction at sea 
level up to some tenths of degrees. Therefore a wind offset 
estimate is computed which provides an angle correction 
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Fig. 7. Pattern-eight flight data comparing the YAE output 4's '° '''^ 
reference value 0r = t^s + {t— '/'w) computed from ship sensors only. 
Except for a drifting offset these curves are in good agreement. This offset 
is compensated in the (/)gr value. For details see Section [vT| Likewise the 
YAE output 9g coincides with the towpoint value 0^. Minor differences 
are due to additional excitation modes of the towing rope and control pod 
system not included in the simple tether assumption. 
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Fig. 8. Flight data comparing sin ipm to the corresponding term of the 
equation of motion (9). We give the unfiltered signal which is a bit noisy 
due to the derivation of ipyn with respect to time. 



value to (p,-. Details of this algorithm go beyond the scope 
of this paper and will be published elsewhere. 

VII. Experimental Validation 

In this section we compare the output angles of the YAE 
to other independent sensor measurements. As already dis- 



cussed in Section VI during tethered dynamic flight the IMU 
and the towpoint measurement values coincide: cpg^. versus 
(ps + {iT — (pw) and dg versus ds- Corresponding measurement 
curves are shown in Fig. IT] The respective curves are in good 
agreement as explained in the figure caption. A verification 
of the ip^ angle output of the navigation setup is shown in 
Fig. Is] The diagram plots the left-hand and right-hand sides 
of (J9k The analogy of the curves shows the validity of this 
equation and thus the capability of the setup in providing a 
reasonable i/'m value. 
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Fig. 9. Estimation algorithm combining YAE and wind referencing. The 
orientation input to the lowpass filter is based on the acceleration vector 
Us and the wind direction Ws which is computed based on the horizontal 
reference ipr and the estimator output angles 9g and i/'g ■ The function blocks 
are to represent the transfer dynamics rather than a concrete algorithm. 



VIII. Discussion and Further Development 

Although the presented YAE shows convincing results 
and has been used successfully for operational flights for 
years, the approach has certain shortcomings in estimating 
the offset rates ujq. The fact that the referencing is with 
respect to the gravity axis only, may lead, in certain specific 
flight situations, to inaccurate offset estimations for sensor 
axes oriented to the Cz-axis. This issue has also been reported 
for a Kalman filter approach on IMU data [11]. In order 
to overcome this situation we suggest a merging of YAE 
and wind referencing into one algorithm, thus processing 
direction of gravity and wind reference simultaneously. 

An approach of a complementary filter processing ac- 
celerometer and magnetometer information is described in 
[14] which could be adopted in the following way: instead 
of the magnetometer we introduce a fictitious sensor value Wg 
emulating a measurement of the horizontal reference wind 
direction. It is obtained by rotating the reference axis e^ into 
the sensor frame by: 



uJs-i?(0r,eg>g)' 



(19) 



Note that 4'g ^nd 9g are filter outputs whereas cj),- is the 
horizontal referencing input. The principal data flow of such 
an algorithm is sketched in Fig. l9] 

We would like to point out here that we consider the 
filtering of the acceleration input as critical: our YAE filter 
is based on the physics that the average of the acceleration 
vectors yields the direction of gravity for a tethered system. 
The common approach first calculates a direction out of the 
acceleration vector dg by normalization and uses as/|as| 
as input value for the filter This is also done for filter 
on the SO(3) described in [14]. As we have to deal with 
accelerations multiple of the gravity acceleration |as| ^ g 
during dynamical flight, it is questionable whether the error 
introduced by neglecting the jag | -amplitude is acceptable. 
Modifications and extensions to these filters may be nec- 
essary. 

First validation tests of the combined algorithms with 
recorded flight data yielded estimated orientation angles 



comparable to the output of the operational YAE and wind 
referencing setup. Meanwhile the estimated sensor biases 
seem to be more realistic than those of the YAE. However a 
critical sensitivity of the bias estimator on filter frequencies 
and flight situations has been observed, which necessitate fur- 
ther research and development efforts on these approaches. 
Conclusively we would like to note that gains and filter 
frequencies of the presented algorithms were determined 
mainly by trial and error on recorded flight data sets. On 
the other hand a thorough Kalman filter design would be 
desirable and may result in some improvements of the esti- 
mator performance. However, such a Kalman theory based 
filter design for a tethered kite system turns out to be non- 
trivial as measurements are less corrupted by sensor noise 
than by disturbances (mainly wind), which can not be easily 
modeled within a statistically based framework. 
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